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Perturbation theory of multi-plane lens effects in terms of mass 
ratios: Approximate expressions of lensed-image positions for two 

lens planes 

Koji IzuMi and Hideki AsADA 
Faculty of Science and Technology, Hirosaki University, Hirosaki 036-8561, Japan 

Continuing work initiated in an earlier publication (Asada, MNRAS 394, (2009) 818), 
we make a systematic attempt to determine, as a function of lens and source parameters, the 
positions of images by multi-plane gravitational lenses. By extending the previous single- 
plane work, we present a method of Taylor-series expansion to solve the multi-plane lens 
equation in terms of mass ratios except for the neighborhood of the caustics. The advantage 
of this method is that it allows a systematic iterative analysis and clarifies the dependence 
on lens and source parameters. In concordance with the multi-plane lensed-image counting 
theorem that the lower bound on the image number is 2^ for TV planes with a single point 
mass on each plane, our iterative results show how 2^ images are realized. Numerical tests 
are done to investigate if the Taylor expansion method is robust. The method with a small 
mass ratio works well for changing a plane separation, whereas it breaks down in the inner 
domain near the caustics. 

§1. Introduction 

Gravitational lensing has become an important subject in modern astronomy 
and cosmology.'^ It has many applications as gravitational telescopes in various 
fields ranging from extra-solar planets to dark matter and dark energy at cosmo- 
logical scales."^ For instance, it is successful in detecting extra-solar planetary sys- 
tems.l3''l3''^''Sl''Sl' Gaudi et have found an analogy of the Sun-Jupiter-Saturn 

system through lensing. Recently gravitational lensing has been used to constrain 
modified gravity at cosmological scale.'^ 

This paper considers the gravitational lensing by point-mass systems on multiple 
planes, where the number of planes is arbitrary. Such a multi-plane treatment is 
important. In microlensing studies, we usually assume a binary lens on a single lens 
plane. In order to discuss its validity, we can consider two lens planes and later 
take the limit that two lens planes merge. In this way, it will become possible to 
estimate the effect caused by a separation between the double lens planes. Another 
importance is for gravitational lensing in cosmology. Clearly, galaxies at different 
redshifts and dark matter inhomogeneities must be described by not a single-plane 
but multi-plane method. 

It has long been a challenging problem to express the image positions as functions 
of lens and source parameters.'^''^ For this purpose, we present a method of Taylor- 
series expansion to solve the multi-plane lens equation in terms of mass ratios by 
extending the previous single-plane work.^^^ In particular, we carefully investigate, 
as a non-trivial task, the denominators of the lens equation with singular points. 

The multi-plane lensed-image counting theorem states that the lower bound on 
the image number is 2^ for planes with a single point mass on each plane (page 
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458 in Fetters, Levine and WambsgansJ^ and references therein). However, the 
counting theorem tells nothing about the image positions. Therefore, it is important 
to discuss how such image positions are realized in an analytical method. 

Under three assumptions of weak gravitational fields, thin lenses and small de- 
flection angles, gravitational lensing is usually described as a mapping from the lens 
plane onto the source plane.l^ Bourassa and KantowskP^'l^^ introduced a complex 
notation to describe gravitational lensing. Their notation was used to describe lenses 
with elliptical or spheroidal symmetry.l^^'l^'l^ 

For point lenses, Witt^ succeeded in recasting the lens equation into a 
single-complex- variable polynomial. This is in an elegant form and thus has been 
often used in investigations of point-mass lenses. The single- variable polynomial due 
to point lenses on a single plane has the degree of N'^ + 1, though the maximum 
number of images is known as 5(A^— l).'23*2iJ'^,)'32j This means that unphysical roots 
are included in the polynomial (for detailed discussions on the disappearance and 
appearance of images near fold and cusp caustics for general lens systems, see also 
Fetters, Levine and Wambsgans^^ and references therein). Following Asada,— ^ we 
consider the lens equation in dual complex variables, so that we can avoid inclusions 
of unphysical roots. 

This paper is organized as follows. In Section 2, the formulation of multi-plane 
lens systems with complex variables is briefly summarized. The lens equation is 
iteratively solved. In section 3, we present iterative solutions for a two-plane case 
and give an algorithm for computing image positions for an arbitrary number of lens 
planes in terms of mass ratios. In section 4, we discuss how lensed-image positions 
are realized in the present method. Section 5 presents numerical tests. Section 6 is 
devoted to the conclusion. 



§2. Basic Formulation 



2.1. Multi-plane lens equation 

We consider lens effects by N point masses, each of which is located at dif- 
ferent angular diameter distances Di {i = 1,2, •••A^) from the observer, where 
Di < D2 < • • • < -Dat. For this case, we prepare N lens planes and assume the 
thin-lens approximation for each lens plane.^^''^^' 

All the deflectors line up in small angles and all are far away from caustics. 
Note that the above lensing setup is idealized. In the real universe, it is rare to find 
a single isolated lensing mass on each plane. We consider that contributions from 
masses at large angles are taken into account in the definition of the angular distance 
on average.'^ In other words, we focus on effects by a local mass distribution 
along the line of sight. 

First of all, angular variables are normalized in the unit of the angular radius of 
the Einstein ring as 



where we put the total mass on the first plane at Di, G denotes the gravitational 
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constant, c means the light speed, Mtot is defined as the total mass YliLi ^'^d 
Di, Ds and Dis denote angular diameter distances between the observer and the 
first mass, between the observer and the source, and between the first mass and the 
source, respectively. 

Recursively one can write down the multi-plane lens equation .1251 ^123 In the 
vectorial notation, the two-plane lens equation is written as 



+ ^2^2 WZTP 



(2.2) 



where (3, 6, l\ and £2 denote the positions of the source, image, first and second 
lens objects, respectively. Here, Vi denotes the mass ratio of each lens object as 
Vi = MiM^l, and we define d2 and 62 as 

It is convenient to use complex variables when algebraic manipulations are done. 
In a formalism based on complex variables, two-dimensional vectors for the source, 
image and lens positions are denoted asw = j3x+ifiy, z = Ox+iOy, and = lix + i^y, 
respectively. Figure [1] shows our notation for the multi-plane lens system. Here, z is 
on the complex plane corresponding to the first lens object that finally deflects light 
rays and thus z means the direction of a lensed image. 

By employing the complex formalism, the two-plane lens equation is rewritten 

as 



w 



\ 



1 — u udo 1 



where the asterisk * means the complex conjugate and we use the identity as 1^1+1^2 = 
1 to delete and i' denotes 1^2- Note that we choose the center of the complex 
coordinate as the first mass position. Then, we have ei = and simply denote 
e = 62, which is the projected relative position of the second mass with respect to 
the first one. The lens equation is non-analytic because it contains not only z but 
also z*. 

2.2. Iterative solutions 

The mass ratio does not exceed the unity by its definition. Therefore, we use a 
simple method of making expansions in terms of the mass ratios. One can delete 1^1 
by noting the identity as Ui = 1. 
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I St Lens Plane 
i Im 



2nd Lens Plane Source Plane 



-f Im 



i Im 




Re 



Fig. 1. Notation: The source and image positions on complex planes are denoted by w and z, 
respectively. Locations of N masses are denoted by for i = 1, ■ ■ ■ , N . Here, we assume the 
thin lens approximation for each deflector. The angular diameter distances among the observer, 
source and each lens object are also defined. 



Formal solutions are expressed in Taylor series as 



oo oo oo 

^= Yl ^s'^^r • • • ^n'hp2){p3)-{pn)^ (2-6) 

P2=0p3=0 pjv=0 



where the coefficients •Z{p2)(p3) - (pjv) independent of any i/j. What we have to do 
is to determine each coefficient -2(p2)(p3)- -(pjv) iteratively. 

At the zeroth order, we have always a single-plane lens equation as the limit 
of 1^1 —> 1 (z^2 = • • • = z/jv — > 0). We have two roots for it. In addition, we have 
other roots for a multi-plane lens equation as seeds for our iterative calculations. An 
algorithm for finding these solutions is explained in next section. 

Note that the above successive approximation cannot work well in the neighbor- 
hood of the caustics, where the mapping becomes singular. Therefore, we focus only 
on regular regions except for the singular domains. 
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§3. Image Positions 

3.1. Two lens planes 

At the zeroth order in u, the two-plane lens equation becomes simply 

w = Z(^o)-^, (3-1) 

^(0) 

which is rewritten as 

^(0)^*0) - 1 = ^^(*o)- (3-2) 
The L.H.S. of the last equation is purely real so that the R.H.S. must be real. Unless 
w = 0, therefore, one can put z^q) = Aw by introducing a certain real number A. By 
substituting 2;(o) = Aw into Eq. (j3-2p . one obtains a quadratic equation for A as 

ww*A^ - WW* A -1 = 0. (3-3) 

This is solved as 




(3-4) 

which gives Z(o) ^ A±w. 

In the particular case of w = 0, Eq. (j3-2p becomes |2;(o)l — 1> which is nothing 
but the Einstein ring. In the following, we assume a general case oi w ^ 0. 

Regarding the denominator of Eq. (|2-5p . we make an expansion in v as 

oo 

zz* - e*z - (1 - i/)52 = Yl ^''/P' (3-5) 

p=0 

where we formally obtain 

/o = ^(o)^(*o) - (^*Ho) - ^2, (3-6) 
/i = Ho)^h + ^{i)^{*o) - + "^s, (3-7) 
/2 = 2(o)^{*2) + ^(i)^(*i) + ^(2)2;(*o) - e*^(2)- (3-8) 

By choosing Z(o) &s A±w, the two-plane lens equation becomes at 0{iy) 

2;(i) + az(\) = 6i, (3-9) 

where we define 

1 
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(3-11) 



The above equation is linear in Z(i-) and thus easily solved as 

bi - abl 
^ ' 1 — aa* 

Next, we consider the two-plane lens equation at 0(z^^). It is written as 

Z(2) +a2;(*2) = 62, (3-12) 

where we define 



(1) 



^(0) 



d2 

fo 



+ — I Zt^i) - Z(o) 



This equation is linear in Z(2) ^^'^ thus easily solved as 

62 - 



Z(2) 



1 — 00* 



(3-13) 



(3-14) 



Let us move to 0(i^^), for which the two-plane lens equation is linearized as 



(3-15) 



where we define 



z, 



7o 



(0) 



^(2) -2(1)-^ +2(0) 



Z(2) + 



7* V 



^(0) 



This equation is easily solved as 



ab% 



2(3) 



1 — aa* 



(3-16) 



(3-17) 



In the above, we have considered a rather general case that the last term in the 
two-plane lens equation (j2-5p is not divergent. Let us investigate the remaining case 
that the denominator of the last term vanishes, which is expressed as 



zz* - (1 - U)52 



€ Z. 



(3-18) 



The left hand side is real and hence the right hand side must be real. Therefore, 
substituting z = ke for a real number k into the above equation leads to a quadratic 
equation for k as 



k'-k-^-^^ = 0. 
ee* 



This is solved to obtain z = z^ 01 z = z^. Here, we define 



z+ 



e + \/e^ + 4z^i<^2e(e*)~^ 
2 



(3-19) 

(3-20) 
(3-21) 
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where they are not expanded in u as an effective renormalization. The reason for 
avoiding an expansion in the denominator is as follows. If an expansion in v were 
done in the denominator, we would see a third-order pole (or higher one) in the 
lens equation. This would lead to more complicated iterations. In order to avoid it, 
therefore, we do not expand the denominator. We exactly treat it. 

It should be noted that the zeros of the denominator are the positions of the lens 
objects in a single-plane lens equation ,'^31' whereas the zeros for the present case are 
not the lens positions but located near the lens positions with a certain correction 
due to a separation between the planes. 

We consider the particular case as z^g) = z±. Then, regarding the denominator 
of Eq. (|2-5p , we make an expansion around z± in v as 

oo 

zz* - e*z - (1 - u)62 = lyPgp, (3-22) 

p=0 

where we formally obtain 

90 = z±z% - e*z± - (1 - v)52 

= 0, (3-23) 

91 = z±zli) + 2(1)4 - (3-24) 

92 = Z±zl2) + ^(l)2(*l) + Z(2)Z± - 6*2(2) • (3-25) 

Here, gp is linear in Zp and z* {p = 1,2, • • • ). Note that a v-tevm appears in Eq. 
()3-23p . because the denominator is exactly treated as a quadratic function. 
At the lowest order in v, the lens equation becomes simply 

w = z± - ^ -—. (3-26) 
This takes the form of 2;(x) + «±2(^-) = 6±i and immediately gives the solution as 

where we define 

«± = (3-28) 

z^- e 

6±i = -. (3-29) 

(4 - e*){w -z± + —) 
z± 

At the next order, the lens equation in the complex-conjugated form is written 

as 

= 2±(i) + ^4(1) + ^ - — ( 2±(i) - z±—) , (3-30) 
^ ' 4 51 V 91/ 
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which is solved as 



Here, we define 

1 



b±2 - a±bl2 .oon 



b±2 



+ 

V / 

(3-32) 

Similarly, we find Z(3). 

Table [7] shows a numerical example of image positions obtained iteratively and 
their convergence. 

tells us an order-of-magnitude estimate of the effect by a separation between 
the two lens planes. Such a depth effect is characterized by ^2, which enters the 
iterative expressions through 2;+ and z_. 

3.2. Three (or more) planes 

The above procedure for two lens planes does not seem to work for an arbitrary 
number of planes, because fifth-order (or higher order) polynomials cannot be solved 
algebraically as shown by Galois.'^ 

By iterative procedures, however, one can construct roots that are nothing but 
image positions, because we have expansion parameters. Let us explain this iterative 
calculation for three-plane lenses as a simple example. One can write down the three- 
plane lens equation. 

First, we neglect 1^3 terms, so that the three-plane lens equation can be reduced to 
the two-plane one. We have already known how to construct four functions denoting 
image positions for the two-plane lens equation. Next, in the similar manner to the 
two-plane case, one can substitute the perturbative image positions into the three- 
plane lens equation. Four positions with the correction at 0(1^3) are thus obtained. 
By using these four linear-order roots, one can find four image positions at 0(z^|). 
In this way, one can recursively obtain higher order roots. 

Other image positions come from the denominator of the last term of the three- 
plane lens equation. The denominator takes the same form as the two-plane lens 
equation but with different coefficients (obtained by a replacement as 5" ^ 3 in 
the subscripts). Hence one can perturbatively construct four other roots. By using 
these four roots as seeds for further iterations, one can construct four roots that can 
perturbatively satisfy the three-plane lens equation. 

Therefore, one can perturbatively construct totally 4 + 4 = 8 image positions. 
Clearly this procedure can be used also for four-plane lens systems. 

First, we ignore 1^4 terms in the four-plane lens equation, so that the equation 
can be reduced to the three-plane lens equation. For = 3, one can find eight 
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image positions as discussed above. Hence, one can iteratively obtain an iterative 
expression of eight image positions in terms of 1/4. Next, let us take a look at the 
denominator of the 1/4 term in the four-plane lens equation. Finding roots of the 
denominator is essentially similar to that of the three-plane lens equation. This can 
be done. More eight roots are thus obtained as seeds for iterative calculations. One 
can perturbatively construct totally 8 + 8 = 16 image positions as functions of the 
source and lens parameters. 

We continue the iterative procedure for five (or more) lens planes step by step, 
so that image positions can be perturbatively obtained as functions of the source 
and lenses. 

Note that two images can merge in the vicinity of the caustics. In this paper, 
we consider only the regular regions, where images cannot merge. Therefore, zeros 
of the denominator of the lens equation are not degenerate but distinct. 

Table I. Example of image positions by the two-plane lens. We choose vi — 9/10, ^2 = 1/10, 
e = 3/2, w = 2, Di/Ds = 2/5, D2/DS = 3/5. Iterative results (denoted as '0th'. '1st', 
'2nd' and '3rd') show a good convergence for the value (denoted as 'Num') that is obtained 
by numerically solving the lens equation. For the same parameter value, a simple ray-tracing 
method gives numerical values (in the row denoted as 'Ray'). 



Images 


1 


2 


3 


4 


0th. 


2.414213 


-0.414213 


1.780776 


-0.280776 


1st. 


2.434312 


-0.390217 


1.731605 


-0.276050 


2nd. 


2.430981 


-0.388713 


1.732327 


-0.275043 


3rd. 


2.431474 


-0.388781 


1.732190 


-0.274861 


Num 


2.431396 


-0.388766 


1.73220 


-0.274833 


Ray 


2.i:!2 


-0.393 


i.732 


-0.271) 



§4. Realizing Images for N Point Masses 

Instead of seeking explicit expressions of image positions, in this section, we 
discuss how to perturbatively realize lensed-image positions for arbitrary N planes. 
A hint has appeared in the previous section. 

For N = 2, the number of the images that are obtained perturbatively is four, 
which equals to 2^^ for N = 2. By induction, we shall show how at least 2-'^ images 
are realized for N lens planes except for the neighborhood of the caustics. 

Let us assume that at least 2^ images are realized for N = p. Note that they are 
not degenerate, since we do not consider the neighborhood of the caustics. What we 
have to do is to show that at least 2*'+^ images are realized for A'" = p -|- 1. 

We consider p + 1 lens planes. First, let us ignore the (p 4- l)-th mass term in 
the lens equation, so that the equation has the same structure as that for N = p. 
By the assumption for N = p, therefore, the reduced equation with neglecting the 
(p + l)-th mass term has at least 2^ roots. 

Next, the (p-l-l)-th mass can be considered a new perturber. In the lens equation, 
the denominator of the fraction with i/p+i has 2^ zeros. Note that it cannot be 
factored because it is a polynomial mixed with z and z*. These zeroth-order roots 
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as seeds lead to iterative image positions with the same number. 

In total, at least 2^ + 2^ = 2^+^ roots are realized, since our iteration method 
does not exclude additional solutions. By induction, we understand how at least 2^ 
images are realized for the multi-plane lens equation for arbitrary N except for the 
neighborhood of the caustics. 

The above method of constructing the image positions means that the image 
number inequality > 2^ is sharp and the lower bound is actually attained. Obstruc- 
tion points at which the backward-traced light ray hits a lens object and hence does 
not reach the source plan^^ play a role in the realization in the sense that some of 
the images are found by investigating the neighborhood of obstruction points. 

Before closing this section, we make two remarks. The first remark is made 
upon a comparison with numerical methods The present method gives analytical 
expressions of image positions (not their value but their functional forms), so that 
there can be two merits: 1) Calculations are faster when we obtain the numerical 
image position. 2) Dependence on the parameters can be made clearer. However, 
it has a disadvantage, because the result seems to take very lengthy expressions. 
The second remark is made on multiple roots. The Taylor method assumes that the 
Taylor series converges, whereas at multiple roots (i.e. merging images) it becomes 
divergent. Therefore, it is unlikely that the Taylor method can fix the problem when 
images merge. See also the next section. 



We perform simple ray-tracing calculations in order to investigate if the Taylor 
expansion method is robust. Table U] shows that numerical results by both methods 
of the Taylor expansion and the ray tracing are in agreement. Note that especially 
the image No. 1 and 3 are in good agreement, though the image No. 2 and 4 have a 
few percent diff'erence. This is because light rays corresponding to images No. 2 and 
4 pass closer to the primary lens compared with No. 1 and 3 and therefore numerical 
errors become relatively large. 

We make also numerical tests for various values of the lensing parameters in 
order to investigate the typical size of the cross section, plane separation and mass 
ratio for which the Taylor expansion breaks down. Figure [2] shows the accuracy 
when the lensing parameters are numerically changed. Here the relative error for 
each image for the chosen parameters is denoted as 



where ZTayior denotes a root obtained by the Taylor expansion method (including the 
third order corrections) and ZNum denotes the root that is obtained by numerically 
solving the lens equation. AMax denotes the largest error among four image positions 
(for two planes) for the chosen parameter values. For Figure [21 we assume the 
parameter values that are the same as those in Table HI namely Di = 0.4, D2 = 0.6 
{D12 = 0.2), v = 0.1 (mass ratio), e = 1.5 (secondary mass position), w = 2 (source 
position) . The Taylor expansion method is robust for a wide range of the parameters 



§5. Numerical tests 



A = 



^Taylor ^NuTn 




^Num 
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as shown by Figure [2j However, it does not work well in some cases, for instance 
when a mass ratio is large (i.e. comparable masses), or the source is located in the 
inner region near the caustics. 

In order to quantify the breakdown of the present method, we choose the thresh- 
old for AMax as 0.01 (one percent). The typical size of the cross section for which 
the Taylor expansion method breaks down is approximately TT X 0.2^ ~ 0.1 (near the 
primary mass direction) and vr x 1^ ~ 3 (around the secondary one), respectively. If 
the threshold is less than one percent, the two cross sections are merged as a single 
one. Note that the angle is normalized by the Einstein radius with the total mass 
located at Di. The typical size of the mass ratio to invalidate the Taylor expansion 
is 0.3. As for the plane separation, the Taylor expansion method seems robust. 



§6. Conclusion 



We made a systematic attempt to determine, as a function of lens and source 
parameters, the positions of images by multi-plane gravitational lenses. We presented 
a method of Taylor-series expansion to solve the multi-plane lens equation in terms 
of mass ratios except for the neighborhood of the caustics. 

In concordance with the multi-plane lensed-image counting theorem that the 
lower bound on the image number is 2^ for planes with a single point mass on 
each plane, our iterative results directly show how 2^ images are realized except for 
the neighborhood of the caustics. 

It is left as a future work to compare the present result with state-of-art numer- 
ical simulations. 
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Fig. 2. Numerical tests of the accuracy of the Taylor expansion method with different parameter 
values. As a reference model for comparisons, we choose the model parameters as ui — 9/10, 
1/2 ~ 1/10, e = 3/2, w = 2, Di/Ds ~ 2/5, D2/DS ~ 3/5, which are the same as those in Table 
1. Top: Mass ratio as 1/ = V2 is changed. Middle: Plane separation as D12 is changed. Note that 
£'12 < 0.6 since Di — 0.4 (normalized by Ds)- Bottom: Source position w is changed along the 
real axis. In actual calculations, smaller parameter steps are adopted to investigate the regions 
near the primary and secondary caustics. The vertical axis denotes the largest relative error 
^Max of the four images. 
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